Correlation functions and momentum distribution of one-dimensional Bose systems 
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The ground-state correlation properties of a one-dimensional Bose system described by the Lieb- 
Liniger Hamiltonian are investigated by using exact quantum Monte Carlo techniques. The pair 
distribution function, static structure factor, one-body density matrix and momentum distribution 
of a homogeneous system are calculated for different values of the gas parameter ranging from the 
Tonks-Girardeau to the mean-field regime. Results for the momentum distribution of a harmonically 
trapped gas in configurations relevant to experiments are also presented. 
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The recent experimental achievements in realizing 
quasi one-dimensional (ID) quantum degenerate Bose 
gases confined in harmonic traps [1-3], have attracted 
a great interest on the ground-state properties of these 
systems. The ID regime is reached in highly anisotropic 
traps, where the axial motion of the atoms is weakly con- 
fined while the radial motion is frozen to zero point os- 
cillations by the tight transverse confinement. In these 
conditions and if the characteristic range of interparticle 
interactions is much smaller than the typical length of the 
transverse confinement, the system can be described by 
the Lieb-Liniger model of ^-interacting ID bosons. For a 
homogeneous system the Hamiltonian has the form 
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where m is the atomic mass and gm = 2Ti 2 /m\a\D\ is 
the coupling constant in terms of the ID effective scatter- 
ing length a id [4]. Recent numerical simulations using 
exact quantum Monte Carlo techniques have explicitly 
shown the cross-over from the three-dimensional (3D) to 
the ID Lieb-Liniger regime for the energetics and struc- 
tural properties of trapped systems [5,6]. The Hamil- 
tonian Hll is expected to provide the correct descrip- 
tion also for the correlation properties of a quasi ID 
system. In this Letter we calculate for the first time, 
using exact quantum Monte Carlo methods, the pair dis- 
tribution function and the one-body density matrix of 
a Lieb-Liniger gas in its ground state. By taking the 
Fourier transform of the one- and two-body correlation 
function we obtain results for the momentum distribution 
and the static structure factor of the system respectively. 
These correlation functions can be accessed in experi- 
ments using ballistic expansion and Bragg spectroscopy 
techniques and can provide a signature of the exotic prop- 
erties of these quantum degenerate ID systems. 

The equation of state and excitation spectrum of a 
homogeneous system of bosons described by the Hamil- 
tonian (1) have been calculated exactly in the thermody- 
namic limit in Rcf. [7] for any value of the dimensionlcss 
gas parameter n|airj|, where n is the particle density. 



If n|ai£)| 3> 1, the system is in the weakly interacting 
regime and mean-field theory can be applied. In this case, 
the energy per particle is given by the Gross-Pitaevskii 
(GP) prediction: E/N — girjn/2. In the opposite limit, 
n|ai£> | <C 1, the system enters the Tonks-Girardeau (TG) 
regime of a gas of impenetrable bosons. In this regime, 
which has been first investigated by Girardeau [8], the 
system acquires fermionic properties in the sense that 
there exists an exact mapping between the wave function 
of the interacting bosons and the wave function of a non- 
interacting Fermi gas. The energy per particle coincides 
in this case with the Fermi energy: E/N = ir 2 h 2 n 2 /(6m). 

We study the following correlation functions: the one- 
body density matrix, which in terms of the ground-state 
wave function ^ (zi, .., zjv) is defined as 



gii z ) = 



N J *g(zi + z, .., z N )^p(zi, ..,zn) dz 2 ..dzN 
n J\^f (z 1 , -,zn)\ 2 dz\..dz N 
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and the pair distribution function 

N(N - 1) / |*o(2:i, .., z N )\ 2 dz 3 ..dz N 
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n 2 J |*o(2i, .., z N )\ 2 dz\..dz 



N 



(3) 



By taking the Fourier transform of these correlation func- 
tions one obtains the momentum distribution n(k) — 
n J dz e lkz gi(z) and the static structure factor S(k) = 
l + njdz e lkz [g 2 (z) - 1]. 

For a system of impenetrable bosons (n|ai£>| <C 1) 
the above correlation functions can be calculated analyt- 
ically. By exploiting the Bose-Fermi mapping one finds 
the following result for the pair distribution function [8] 



g 2 {z) = l-jl(im\z\) , 



(4) 



where jo is a spherical Bessel function. The correspond- 
ing static structure factor is given by 
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The one-body density matrix of a TG gas has also been 
calculated by performing analytic expansions for short 
and long distances (see Ref. [9] and references therein). 
The long-range behavior is given by g\{z) oc 1/ \/\z\ [10], 
corresponding to an infrared divergence in the momen- 
tum distribution n{k) cx l/-y/|fc|. 

Beyond the TG regime full closed-form expressions of 
the correlation functions are not known. The long-range 
asymptotics of the one-body correlation function can be 
obtained from the hydrodynamic theory of low-energy 
excitations [11,12]. One finds the following power-law 
decay 



gx{z) oc l/\z\ 



(6) 



where a = 7nc/{2irhn) is fixed by the density n and by 
the speed of sound c. The above result holds for \z\ 3> £ 
where £ = Tij{\f2mc) is the healing length of the system. 
The velocity of sound can be obtained from the equation 
of state through the compressibility of the system mc 2 = 
nd^i/dn, where /i = dE/dN is the chemical potential. In 
the TG regime (n|ai£>| -C 1) one finds mc = nhn and 
thus a = 1/2 as anticipated above. In the opposite GP 
regime (n|ai£>| ^> 1), the result is a — l/(7r^/2n|ai£)|) 
which decreases as n|air)| increases. Of course, result 
(6) excludes the existence of Bose-Einstein condensation 
even in the ground-state [13]. The infrared behavior of 
the momentum distribution follows immediately from Eq. 
(6) 



n(ife) oc l/\k\^ 
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holding for \k\ <C l/£. Furthermore, the hydrodynamic 
theory allows one to calculate the static structure factor 
in the long-wavelength regime |fc| <C l/£. One finds the 
well-known result 



S(k) -» Ti\k\/{2mc) , 



(8) 



characteristic of phonon excitations. Recently the short 
range behavior of the one- and two-body correlation func- 
tions has also been investigated. The value at z = of 
the pair correlation function can be obtained from the 
equation of state through the Hcllmann-Fcynman theo- 
rem [14]: g 2 (z = 0) = -(m\a 1D \/h 2 )n 2 d[(E/N)/n 2 }/dn. 
The result approaches zero in the TG regime and unity 
in the GP regime. Furthermore, the first few terms of the 
short-range series expansion of the one-body correlation 
function g\{z) = 1 + X^i c i\ nz \ l nave been calculated in 
Ref. [9] again exploiting the knowledge of the equation 
of state. 

In this Letter we provide a complete calculation of 
the spatial dependence of the one- and two-body corre- 
lation functions for values of the gas parameter ranging 
from the TG to the deep GP regime. The calculation 
is carried out using the diffusion Monte Carlo (DMC) 
method, which allows one to solve exactly apart from 



statistical uncertainty, the many-body Schrodinger equa- 
tion for the ground state of a Bose system [15]. We con- 
sider a system of N particles described by the Hamilto- 
nian (1), with periodic boundary conditions. The sys- 
tem has size L and linear density n = N/L. Impor- 
tance sampling is used through the trial wave function 
if) T {zi,..,z N ) = n»<j /(%')> where f(z) is a two-body 
term which we choose as 



/(*) 



Acos[k(\z\ - B)\ 
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(\z\>Z). 
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The z = boundary condition .f(0+) - /' (0~) = 
2/(0)/|ai£>|, which accounts for the ((-function po- 
tential, fixes the parameter k through the condition 
fc|ai£>| tan(A;.B) = 1. The remaining parameters A, B 
and [3 are fixed by requiring that the function f(z), its 
derivative f'(z) and the local energy —2f"(z)/f(z) be 
continous at the matching point z = Z. The value of 
the matching point Z is a variational parameter that 
we optimize using variational Monte Carlo (VMC). For 
\z\ < Z the Jastrow term corresponds to the exact solu- 
tion of the two-body problem and is expected to provide 
a correct description of short-range correlations. Long- 
range correlations arising from phonon excitations are 
instead accounted for by the functional dependence of 
f(z) for \z\ > Z [11]. Notice that at \z\ = L/2, which 
sets the limits of our simulation box, / = 1 and par- 
ticles are uncorrelated. The TG wave function *q g = 
Y\ i< j | sin[7r(zi — Zj)/L]\ is obtained as a special case of 
our trial function ipT for Z = B = L/2 and kL = it. 
The optimized Jastrow function introduced above gives 
for all values of n|cii£)| variational results for the energy 
which are almost identical to the exact DMC ones. The 
level of accuracy of the trial wave function is particularly 
important for the calculation of gi(z). In fact, the pair 
distribution function g^{z) is calculated using the method 
of "pure" estimators, unbiased by the choice of the trial 
function [16]. Due to the non-local property of the one- 
body density matrix, the function gi(z) can instead be 
obtained only through the extrapolation technique. For 
an operator A which does not commute with the Hamil- 
tonian, the output of the DMC method is the mixed esti- 
mator ('JoI^IV't) and one can approximate the "pure" es- 
timator by (* |^4|*o) = 2(* |-4|V't) - (ipT\A\ip T ) , where 
(iPt\A\iPt) is the variational estimator. Of course, this 
procedure is very accurate only if ip? — ^o- We find that 
DMC and VMC give results for gi (z) which are very close 
and we believe that the extrapolation technique is in this 
case exact. We also find that gi(z) is strongly affected 
by finite size effects particularly in the GP regime, where 
simulations of quite large systems are needed to obtain 
stable results. In order to control these effects we per- 
form simulations with number of particles ranging from 
N = 100 to N = 500. 

In Fig.l we show the results of the pair distribution 
function for values of n|aiu| ranging from the TG to the 
GP regime (the corresponding values of the energy per 
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particle are shown in the inset) . For the smallest value of 
the gas parameter (n|an)| = 10- 3 ) we find that g^ (z) ex- 
actly coincides with the TG result (4). We also find that 
the values of 32(0) exactly reproduce the results obtained 
from the equation of state. We notice that for small val- 
ues of the gas parameter correlations die out very quickly 
and g2{z) — 1 after few interparticle distances. For larger 
values of n|ai£>| correlations extend over much larger dis- 
tances. This result is also visible from the behavior of the 
static structure factor (shown in Fig. 2) which in the GP 
regime is dominated by long-wavelength fluctuations. 
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FIG. 1. Pair distribution function for different values of 
the gas parameter. In ascending order of the value at z = 0: 
fi|ai£>| = 10 , 0.3, 1, 30, 10 . Arrows indicate the value of 
52(0) as obtained from the equation of state. Inset: energy per 
particle for the same values of n|oi£)| (symbols), Lieb-Liniger 
equation of state (solid line), GP limit (dashed line), TG limit 
(dotted line). Energies are in units of h 2 / \2ma\ D ). 

In Fig. 3 we present results for the one-body density 
matrix. We see that for distances considerably larger 
than the healing length £, g\{z) is well reproduced by 
the asymptotic behavior (6) for which we calculate the 
proportionality coefficient through a best fit to the DMC 
results. The deviations from the power law decay at the 
largest distances (z ~ L/2) are due to finite size effects. 
In Fig. 4 we show the corresponding results for the mo- 
mentum distribution. This is obtained from the Fourier 
transform of the calculated one-body correlation func- 
tion at short distances and of the fits to the power-law 
decay (6) at large distances. We notice that the infrared 
asymptotic behavior (7) is recovered for values of k con- 
siderably smaller than the inverse healing length l/£. At 
large k the numerical noise of our results is too large to 
extract evidences of the 1/fc 4 behavior predicted in [9] . 

Let us now turn to harmonically trapped systems. In 
the ID regime the system is described by the Hamiltonian 
HhL + Y^i=i mLJ z z i/2 w ith an effective coupling constant 
given by gin = 2h 2 a/(ma 2 L ), where a is the 3D s-wave 



scattering length and a± = y/h/(mu>±) is the length 
fixed by the tight transverse confinement [4,6]. Conse- 
quently, the value of the effective ID scattering length 
is given by the ratio |oid| = a\/a. The relevant pa- 
rameters are: the ratio a/a±, the anisotropy parameter 
A = uj z /uj± and the number of trapped particles N. 




FIG. 2. Static structure factor for the same values of n|oiu | 
as in Fig. 1 (solid lines). The dashed lines are the correspond- 
ing long- wavelength asymptotics from Eq. (8). 
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FIG. 3. One-body density matrix for the same values of 
n|oiu| as in Fig. 1 (solid lines). The dashed lines are the 
fits to the corresponding long-wavelength asymptotics from 
Eq. (6). The arrows indicate the value of £n: the leftmost 
corresponds to n|oio| = 10~ 3 , the rightmost to n|oio| = 10 3 . 

The DMC simulation is carried out using the follow- 
ing trial wave function for the importance sampling: 
i> T (zi,..,z N ) = Yli9(zi)Yli<j f( z ij), where 9{z) = 
exp(— a z z 2 ) is a one-body term which accounts for the 
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confinement in the longitudinal direction and is fixed by 
the variational parameter a z . In the two-body term we 
neglect long-range correlations, which are irrelevant in 
the trapped case, and we choose f(z) as given by (9) for 
\z\ < Z and f(z) = 1 for \z\ > Z. We consider the follow- 
ing configurations: a/a± = 0.2, A = 1CP 3 and number of 
particles N =5, 20, 100. In Rcf. [6] we have proven that 
in these conditions the ground-state energy and structure 
of the cloud is correctly described by the Licb-Linigcr 
equation of state in local density approximation. The 
momentum distribution of a trapped system is obtained 
as: n(k) = J dZdz' n{Z)g 1 (Z+z'/2, Z-z'/2)e ikz \ where 
n(z) is the density profile and n[(z + z')/2]gi(z, z') = 
N J .., z N )^ (z', .., z N )dz 2 ..dz N / J |*o(2i, .., zat)| 2 

dz\..dzN- 
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FIG. 4. Momentum distribution for the same values of 
n|ai£>| as in Fig. 1 (solid lines). The dashed lines correspond 
to the infrared behavior of Eq. (7). The arrows indicate the 
value of l/(£n): the rightmost corresponds to n|air)| = 1CF 3 , 
the leftmost to n|anj = 10 3 . 

In Fig. 5 we show the results for the momentum distri- 
bution. Note the peaks becoming sharper with increasing 
N. In the case of N — 5 and N — 20, we find that the 
rounding off of n(k) at k <~ 1/R Z , where R z is the size of 
the cloud in the axial direction, washes out completely 
the divergent behavior at small momenta. For the largest 
system with N = 100 the value of the gas parameter is 
no|aiD| — 1-1, where no is the density in the center of the 
trap, and the corresponding value of the exponent a is 
a ~ 0.19. In the region of wave vectors 1/R Z < k < l/£, 
where the healing length £ is estimated in the center of 
the trap, we find some evidence of the infrared behav- 
ior (7) (sec inset of Fig. 5). In order to see a cleaner 
signature of this effect one should consider much larger 
systems. 

In conclusion, we have carried out a complete study 
of the one- and two-body correlation functions of the ID 
Lieb-Linigcr model at T = 0. We have also investigated 



how the presence of harmonic trapping modifies the be- 
havior of the momentum distribution in configurations 
relevant to experiments. 
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FIG. 5. Momentum distribution of a trapped system. In- 
set: momentum distribution for N = 100 (solid line) on a 
log-log scale. The dashed line is a fit to l/fc 1- " with a = 0.19. 
The momentum distribution is in units of a z = yjh/{mu z ). 
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